NPS-MA-9 7-001 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


On  Satellite  Umbra/Penumbra  Entry 
and  Exit  Positions 

by 

Beny  Neta 

David  Vallado,  LtCol 

January  1997 


o 


Approved  for  public  release;  distribution  is  unlimited. 

Prepared  for:  Office  of  Naval  Research 

Monterey,  CA  93943 

and 

Phillips  Laboratory 
Kirtland  AFB,  NM87117 


mtC  WAUT*  jttiHMSCTBD  1 


NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


Rear  Admiral  M.  J.  Evans  Richard  Elster 

Superintendent  Provost 

This  report  was  prepared  in  conjunction  with  research  conducted  for  Office  of  Naval  Research  and  funded  by 
Phillips  Laboratory. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 

This  report  was  prepared  by: 


Beny  Neta 

Professor  of  Mathematics 


Lt  Col  (sel)  David  Vallado 
Orbital  Analyst 


Reviewed  by: 

JV-  /$• 

WALTER  M.  WOODS 
Chairman 


Associate  Provost  and  Dean  of  Research 


REPORT  DOCUMENTATION  PAGE 

Form  approved 

OMB  No  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  tbe  time  for  reviewing  instructions,  searching  existing  data  sources 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  information  Operations  and  Reports,  1215  Jefferson 

DaYK  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington  DC  20503 

1.  AGENCY  USE  ONLY  (Leave  blank) 

2.  REPORT  DATE 

January  1997 

3.  REPORT  TYPE  AND  DATES  COVERED 

Technical  Report:  1  July  1996-31  December  1996 

4.  TITLE  AND  SUBTITLE 

On  Satellite  Umbra/Penumbra  Entry  and  Exit  Positions 

5.  FUNDING 

FMBD-96-560-NPS 

6.  AUTHOR(S) 

Beny  Neta  and  David  Vailado 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESSES) 

Naval  Postgraduate  School 

Monterey,  CA  93343-5000 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

NPS-MA-97-001 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Office  of  Naval  Research  and  Phillips  Laboratory 

Monterey,  CA  93943  Kirtland  AFB,  NM  87117 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

The  views  expressed  in  this  report  are  those  of  the  author  and  do  not  reflect  the  official  policy  or  position  of  the  Department  of  Defense 
or  the  United  States  Government. 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words.) 

The  problem  of  computing  Earth  satellite  entry  and  exit  positions  through  the  Earth’s  umbra  and  penumbra,  for  satellites  in  elliptical  orbits,  is 
solved  without  the  use  of  a  quartic  equation.  A  condition  for  existence  of  a  solution  is  given.  This  problem  is  related  to  perturbation  force  resulting 
from  solar  radiation  pressure. 

14.  SUBJECT  TERMS 

Umbra/penumbra,  Halley’s  method,  Newton’s  methods 

15.  NUMBER  OF 

PAGES 

24  pages 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION  1  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION 

OF  REPORT  1  OF  THIS  PAGE  OF  ABSTRACT 

Unclassified  |  Unclassified  Unclassified 

20.  LIMITATION  OF 
ABSTRACT 

SAR 

NSN  7540-01-280-5800  Standard  Form  298  (Rev.  2-89) 


Prescribed  by  ANSI  Std  239-18 


On  Satellite  Umbra/Penumbra  Entry 
and  Exit  Positions 

Beny  Net  a* 

Naval  Postgraduate  School 
Department  of  Mathematics 
Code  MA/Nd 
Monterey,  CA  93943 

Lt  Col  (sel)  David  Vallado 
Orbital  Analyst 
Phillips  Laboratory 
Kirtland  AFB,  NM  87117 

December  19,  1996 


‘Author  to  whom  all  correspondence  should  be  addressed. 


1 


Abstract 


The  problem  of  computing  Earth  satellite  entry  and  exit  positions  through 
the  Earth’s  umbra  and  penumbra,  for  satellites  in  elliptical  orbits,  is  solved 
without  the  use  of  a  quartic  equation.  A  condition  for  existence  of  a  solution 
is  given.  This  problem  is  related  to  perturbation  force  resulting  from  solar 
radiation  pressure. 

Keywords:  umbra/penumbra,  Halley’s  method,  Newton’s  method 


1  Introduction 

The  problem  of  computing  Earth  satellite  (in  elliptical  orbits)  entry  and  exit 
positions  through  the  Earth’s  umbra  and  penumbra  is  a  problem  dating  from 
the  earliest  days  of  the  space  age,  but  it  is  still  of  the  utmost  importance  to 
many  space  projects  for  thermal  and  power  considerations  (Mullins,  1991). 
It’s  also  important  for  optical  tracking  of  a  satellite.  To  a  lesser  extent,  the 
satellite  external  torque  history  and  the  sensor  systems  are  influenced  by  the 
time  the  satellite  spends  in  the  Earth’s  shadow. 


Figure  1:  Earth  umbra  and  penumbra 
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The  umbra  is  the  conical  total  shadow  projected  from  the  Earth  on  the 
side  opposite  the  Sun.  In  this  region,  the  intensity  of  the  solar  radiation 
is  zero.  The  penumbra  is  the  partial  shadow  between  the  umbra  and  the 
full-light  region  (see  figure  1).  In  the  penumbra,  the  light  of  the  Sun  is 
only  partially  cut  off  by  the  Earth,  and  the  intensity  is  between  0  and  1. 
All  textbooks  discussing  the  problem  (e.g.  Geyling  and  Westerman,  1971, 
and  Escobal,  1985)  even  the  recent  work  by  Mullins  (1991),  suggest  the 
use  of  a  quartic  equation  analytic  solution.  Because  the  quartic  is  a  result 
of  squaring  the  equation  of  interest,  one  must  check  all  four  solutions  and 
discard  the  spurious  ones.  In  this  paper,  we  examine  solving  the  original 
equation  numerically.  We  will  give  a  condition  for  the  existence  of  a  solution, 
discuss  the  initial  guess  for  the  iterative  scheme,  and  compare  the  complexity 
of  the  two  schemes  (ours  versus  the  analytic  solution  of  the  quartic). 

The  shadow  problem  has  been  solved  in  the  past  by  assuming  a  cylindri¬ 
cal  shadow  behind  the  Earth  (Geyling  and  Westerman,  1971),  or  a  conical 
shadow  which  is  more  realistic  (Fixler,  1964,  and  Mullins,  1991).  The  nu¬ 
merical  solution  will  be  discussed  for  each  case. 

2  Problem  Formulation 

In  this  section,  we  formulate  the  problem  using  both  cylindrical  and  conical 
shadow  geometry.  We’ll  see  that  the  solution  method  is  different  in  the  two 
cases. 


2.1  Cylindrical  Shadow 

In  this  case  the  orbital  geometry  is  given  in  figure  2  (Escobal,  1985,  p.  157, 
or  Vallado,  1996,  p.521). 

The  analysis  given  by  Escobal  (1985)  and  Vallado  (1996)  show  that  the 
true  anomaly,  v,  at  entrance  and  exit  into  the  shadow  satisfies  the  following 
equation: 

1  +  e  cos  i/)2  +  p2(/?i  cos  v  +  /?2  sin  u)2  —  p2  =  0  (1) 

where  R$  is  the  radius  of  the  Earth  (~  6378.136  km),  Re  is  the  Sun’s 
position  vector  (~  696000  km),  p  is  semi-parameter,  and  e  is  the  eccentricity. 
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Figure  2:  Cylindrical  Shadow 


The  remaining  classical  orbital  elements  are  inclination,  i,  longitude  of  the 
ascending  node,  0,  and  the  argument  of  perigee,  ui.  The  parameters  0i  and 
02  are  given  by 


•  P 
Rq 


R ©  •  Q 

Rq 


The  unit  vectors  P  and  Q  are  defined  by 


cos  u  cos  Cl  —  sin  a;  sin  Cl  cos  i 
cos  oj  sin  Cl  +  sin  w  cos  Cl  cos  i 
sin  u)  sin  i 


—  sin  u>  cos  0  —  cos  u>  sin  f2  cos  i  * 

—  sin  u>  sin  Cl  +  cos  u  cos  Cl  cos  i 
cos  lo  sin  i 

For  circular  orbits  and  if  i  —  0, 7r,  P  should  be  redefined  in  a  convenient 
manner  (see  Escobal,  1985). 

2.2  Conical  Shadow 

In  this  case,  one  must  distinguish  between  umbra  (full  shadow)  and  penumbra 
(partial  shadow)  regions.  In  the  umbra  case,  we  must  solve  a  system  of  two 
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nonlinear  equations  (see  Mullins,  1991).  The  first  equation  models  the  surface 
of  the  shadow  cone 

^sh’  ysh’  zsh)  =  y2s  h  +  4  -  (d  -  xsh)2  tan2  <*  = 0  (2) 

where  d  is  the  distance  from  center  of  the  Earth  to  apex  of  shadow  cone 
(~  1.3836- 106  km),  and  a  is  half  angle  of  that  cone  (~  .26412°).  The  second 
equation  describes  the  orbit 

G(x0,y0)  =  )2+  (y)2  -  !  =0  (3) 

where  b  =  a\/l  —  e2.  Because  the  two  equations  are  not  in  the  same  coordi¬ 
nate  system,  we  take  and  rotate  it  to  get  r0.  The  transformation  is  given 

by 


r0  =  ROT3(u)  ROTl{i)  R0T3{{1)  ROT l(-e)  ROT3{ir  -  L )  fsh 

where  e  is  the  mean  obliquity  of  the  ecliptic  (~  23.5°  ),  L  is  the  ecliptic 
longitude  of  the  Sun,  and  ROTl((f> ),  ROT3(4>)  are  rotations  about  the  x,  z 
axis  (respectively)  by  <f>.  If  we  denote  the  transformation  matrix  by  A,  then 

ssh  =  anxo  +  fl2iS/o  (4) 

2/sh  =  a12%0  +  <*222/0  (5) 

zsh  =  ai3^o  +  <*232/o-  (6) 

Notice  that  z0  is  zero  at  the  intersection  of  the  two  equations  (2)-(3).  Because 
only  solutions  with  >  0  are  acceptable(see  figure  2),  we  must  satisfy 

duXo  +  <*2i2/o  >  0.  (7) 

Substituting  (4)-(6)  into  (2),  we  get  the  following  equation 

T\(£o,  Vo)  =  c*o£(j  +  aiVo  +  2qt2£o2/o  +  013X0  +  £*42/0  —  d2  tan2  a  =  0.  (8) 

where 

ao  =  a22  -f  a23  —  a2a  tan2  a 

0:1  —  0*22  "f*  <*23  ~ ~  <*21  0" 
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A 

a 2  =  Ol2®22  +  013^23  —  01X^21  tan  <7 
a3  —  2an0?tan2  <r 
or4  =  2a2idtan2  cr 

This  equation  should  be  solved  with  (3)  and  (7). 

Mullins  (1991)  suggests  solving  (8)  subject  to  (3)  and  (7),  using  a  quartic 
equation  for  x  and  then  checking  each  of  the  four  solutions  with  solutions  of  a 
quadratic  equation  for  y  as  a  function  of  x.  Mullins  admits:  “The  coefficients 
(of  the  quartic)  are  messy  functions  of  the  angles  shown  ...”.  In  section  7, 
we  show  a  better  way  to  solve  the  problem  without  going  through  a  quartic 
equation  and  thus  without  computing  these  “messy  coefficients.” 

In  the  penumbra  case,  Mullins  (1991)  shows  that  (2)  becomes 

F(xsh>  S'sh’  zsh)  =  ylh  +  zlh  ~  (d>  +  xsh)2  tan2  o'  =  0 

where  d!  is  the  distance  from  the  center  of  the  Earth  to  the  apex  of  the  cone 
between  the  Sun  and  the  Earth  (~  1.35849  •  106  km),  and  a'  is  half  angle 
of  that  cone  (~  0.26901°).  This  leads  to  an  equation  similar  to  (8)  to  solve. 
The  idea  presented  in  section  7  will  be  used  here  too. 

3  Complexity  of  Quartic  Solution 

The  problem  (for  cylindrical  shadow)  can  be  solved  using  the  quartic  equation 

A0  cos4  v  +  A\  cos3  v  +  A%  cos2  v  +  A3  cos  u  -1-  A4  =  0  (9) 

analytically  and  then  rejecting  the  spurious  roots  based  on  the  following 
conditions:  The  physical  solution,  should  satisfy  the  original  equation  and 

/?i  cos  v  +  /?2  sin  v  <  0. 

The  coefficients  of  the  quartics  are  given  by: 
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2(A!-A2)(l-/J22)-4/3f/322 


If  the  work  is  done  economically,  one  finds  that  the  number  of  multipli¬ 
cations  and  divisions  required  to  compute  the  coefficients  of  the  quartic  is 
38.  To  find  the  solution  of  the  quartic  requires  64  multiplications/division,  5 
square  roots,  4  cube  roots,  1  arccosine  and  3  cosine  evaluations.  The  cosine 
and  arccosine  evaluations  are  required  only  if  the  discriminant  is  negative, 
see  Abramowitz  and  Stegun  (1965). 

4  Numerical  Solution  for  Cylindrical  Shadow 

To  solve  the  shadow  equation  (1)  numerically,  we  can  use  either  bracketing 
or  fixed-point  type  methods.  In  the  following,  we  describe  only  Newton’s 
and  Halley’s  methods  which  are  of  fixed-point  type.  It  is  first  suggested  to 
check  the  existence  of  a  solution.  First,  rewrite  (1)  as: 

f(x)  =  Ax 2  +  Bx  +  CxV  1  —  x2  +  D  =  0  (10) 

where  x  =  cos  v.  In  order  to  have  a  solution,  we  must  have 

/(—!)/(!)  <  0.  (11) 
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Clearly  equality  means  that  cos  v  is  ±1.  The  strict  inequality  in  (11)  is 
equivalent  to 


1  - 


Re  V 

a(l  +  e)J 


>  ft  >  1  - 


2 


Note  that  there  is  no  condition  on  /?2- 


4.1  Newton’s  Method 

To  solve  a  nonlinear  equation  f(x)  =  0  via  Newton’s  method,  we  require  an 
initial  guess  xo-  Then  an  iterative  procedure  can  be  followed  to  construct  a 
sequence  of  estimates  xn,  by 

_  „  f(xn)  „  _  n  i 

Xn+1  —  %n  r,(  u  n  —  0,1,... 

The  iterative  process  converges  if  either 

|/(i„)|  <  Tolj 


or 

i 


|®n+l  %n  ^  T  olx 


for  given  tolerances.  In  either  case  we  take  as  the  root.  The  convergence 
rate  is  quadratic.  If  the  iterative  process  doesn’t  converge  in  a  certain  number 
of  iterations,  we  stop.  In  this  case  we  suggest  bracketing  methods.  Newton’s 
method  will  diverge  if  we  hit  a  point  where  f'(x)  is  very  small. 


4.2  Halley’s  Method 


Halley’s  method  converges  faster  (third  order  compare  to  second  order  for 
Newton).  The  iterative  process  is 


fM 

/WfW’ 

2  f'{xn) 


n  =  0, 1, . . . 
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4.3  Bracketing  Methods 


In  general,  bracketing  methods  are  slower,  but  they  are  safer,  in  the  sense 
that  convergence  is  guaranteed.  For  example,  the  bisection  method  starts 
with  an  initial  interval  containing  the  root,  [cto?^o]-  The  process  halves  the 

interval  at  every  step.  After  n  iterations,  the  length  of  the  interval  containing 
i  .  6o  —  a0 

the  root  is  —  .  Therefore,  the  number  of  iterations  required  depends  on 

the  length  of  the  initial  interval  and  the  tolerance. 

This  simplistic  method  can  be  modified  by  using  Regula  Falsi  (solving  a 
linear  equation  at  every  step)  or  modified  Regula  Falsi  (which  is  useful  when 
the  curvature  of  /  is  large  enough.) 

For  example,  we  have  solved  (10)  with  A  =  1,  B  =  -2,  C  =  1,  and 
D  =  1.  Newton’s  method  required  5  iterations  for  convergence  to  10~10, 
Halley’s  method  required  4  iterations  and  the  bisection  methods  used  31  it¬ 
erations.  If  we  require  a  more  realistic  accuracy,  let  say  10-6,  then  Halley’s 
method  requires  3  iterations,  Newton’s  requires  4  iterations  and  the  brack¬ 
eting  methods  uses  19. 


5  Initial  Guess 

Because  the  problem  is  to  solve  for  cos  i/,  we  know  that  the  solution,  if  it  ex¬ 
ists,  must  lie  in  the  interval  [—1,1].  For  bracketing  methods  we  suggest  using 
this  interval,  and  for  Newton’s  and  Halley’s  method,  we  take  the  midpoint 
of  the  interval,  i.e.  xq  —  0. 

For  subsequent  crossings  through  the  shadow,  we  can  take  x0  to  be  the 
previous  solution. 


6  Complexity  of  Numerical  Solution 

All  iterative  procedures  require  function  evaluations,  and  some  will  require 
the  evaluation  of  the  first  and  maybe  second  derivative.  The  evaluation  of 
the  function  requires  4  multiplications/divisions  (using  nested  multiplication) 
and  1  square  root.  The  evaluation  of  the  first  derivative  is  accomplished  by  7 
multiplications/divisions  and  1  square  root.  The  second  derivative  requires 
8  multiplications/ divisions  and  1  square  root.  For  one  iteration  of  Halley’s 
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method  we  need  23  multiplications/divisions  and  3  square  roots.  For  one 
iteration  of  Newton’s  method  we  need  12  multiplications/divisions  and  2 
square  roots.  For  the  bisection  method  we  need  5  multiplications/divisions 
and  1  square  root.  If  we  multiply  the  number  of  iterations  by  the  cost  per 
iteration  we  find  that  Newton’s  method  is  the  cheapest  with  48  multipli¬ 
cations/divisions  and  8  square  roots,  then  Halley’s  method  with  69  mul¬ 
tiplications/divisions  and  9  square  roots,  then  bisection  with  95  multipli¬ 
cations/divisions  and  5  square  roots.  In  comparison,  Newton’s  method  is 
cheaper  than  solving  the  quartic  and  it  doesn’t  require  checking  for  spurious 
roots.  Even  Halley’s  method  is  competitive  with  the  analytic  solution  of  the 
quartic.  We  summarize  the  results  in  a  table. 


operation 

mult./div. 

sq.  root 

cubic  root 

trig.  func. 

number  iter. 

Newton 

48 

8 

0 

0 

4 

Halley 

69 

9 

0 

0 

3 

Bisection 

95 

5 

0 

0 

19 

Quartic 

102 

5 

4 

2* 

0 

Table  1:  Operation  count 


T  Numerical  Solution  for  Conical  Shadow 

In  this  section,  we  describe  a  numerical  method  to  solve  (8)  and  (3)  subject 
to  (7).  We  suggest  guessing  an  initial  approximation  x0  and  use  (3)  to  get 
the  corresponding  y0 

y0  =  ±bV  1  -  e2.  (12) 

Because  (3)  is  quadratic,  we  offer  here  the  correct  sign  to  satisfy  (7).  Note 
that  (7)  describes  a  half  plane  whose  boundary  is  a  line  in  figures  3  and  4. 

(13) 

*This  doesn’t  include  checking  for  spurious  roots. 


an 

2/0  = - X0. 

021 
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Figure  3:  an  and  021  have  the  same  sign 


Figure  4:  an  and  <2  21  have  opposite  signs 
Therefore  the  sign  of  the  radical  in  (12)  is  the  same  as  the  sign  of  a2\. 
We  now  rewrite  (8)  as 

F^x,  y )  =  Ar2  +  By{x )2  +  Cxy(x)  +  Dx  +  Ey{x )  +  F 
with  (using  (3)) 

y(x)  =  ±y/l  —  e2^/a2  —  (x  +  ae)2 
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For  Newton’s  method,  we  need  F[  and  y'  which  are  given  by 


F[{x,  y )  =  2Ax  +  2  By(x)y\x)  +  Cy(x)  +  Cxy'(x)  +  D  +  Ey\x) 


y 


and 

,  (1  -  e2)(x  +  ae ) 

y\x)  =  =F~ - - - 1 

Now  the  iterative  procedure  is  as  follows 

Fi(xn,yn) 


xn+l  xn 


(#n>  J/n) 


n  =  1,2, .  • . 


2/n+i  =  ±\/l  -  e2^J a2  -  (xn+i  +  ae)2 

Remember  to  choose  the  sign  appropriately. 


8  Conclusions 

In  this  paper,  we  suggest  the  use  of  iterative  techniques  to  compute  the 
entry  and  exit  positions  through  the  Earth’s  umbra  and  penumbra.  We  also 
show  how  to  choose  the  initial  guess  for  the  first  and  subsequent  crossings. 
Several  iterative  methods  for  the  solution  of  the  problem  are  compared  to  the 
currently  used  one.  Newton’s  method  converges  fast  especially  at  subsequent 
crossings,  because  the  initial  guess  is  close  enough. 
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Appendix  A 


This  section  provides  two  FORTRAN  subroutines  using  Newton’s  method 
to  solve  both  the  cylndrical  and  cone  shadows. 


c 

c  Newton's  algorithm  for  the  cylindrical  shadow 

c 

real*8  xn,fn,xnl,fpn,tol,tolf 
real*8  a0,al,a2,a3 
integer  indx,mxindx 
c 

c  max  number  of  iterations 
c 

mxindx=15 

c 

c  convergence  tolerance  on  consecutive  iterates 
c 

c  tol=l.e-18 
tol=l . e-6 
c 

c  convergence  tolerance  on  closeness  of  function  to  zero 

c 

c  tolf=l.e-18 

tolf=l .e-6 
aO=l . 
al=-2. 
a2=-l . 
a3=l . 
c 

c  initial  guess 

c 

xn=0 . 
indx=0 
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continue 

call  f (a0,al,a2,a3,xn,fn,fpn) 

is  function  close  to  zero  at  xn? 

if (dabs(fn) .le.tolf)  go  to  20 

compute  next  iterate 

xnl=xn-fn/fpn 
znl=dsqrt (1 . d0-xnl*xnl) 
indx=indx+l 

print  *,indx,xn,xnl,znl,fn 

check  for  closeness  of  iterates 

if (dabs(xnl-xn) .le.tol)  go  to  10 
xn=xnl 

check  if  max  number  of  iterates  exceeded 

if (indx.ge.mxindx)  go  to  30 
go  to  5 

print  convergence  iterates  close’ ,xn,xnl 
go  to  40 

print  *, ’  convergence  function  close  to  zero’, xn,fn 
go  to  40 

print  no  convergence  -  max  number  of  iterates’ 

stop 

end 

subroutine  f (a0,al,a2,a3,x,y,yp) 
real*8  x,y,yp 
real*8  a0,al,a2,a3 

evaluate  the  function 


z=dsqrt ( 1 . d0-x*x) 


y=aO*x**2+al*x+a2*x*z+a3 


c 

c  evaluates  the  first  derivative 
c 

yp=2 . *a0*x+al+a2*z-a2*x/z*x 

return 

end 


c 

c  Newton's  algorithm  for  the  cone  shadow 

c 

real*8  xn , yn , f n , xnl ,ynl , ypn , f pn , t ol , tolf 
real*8  all,al2,al3,a21,a22,a23 
real*8  a,b,c,d,e,f 

real*8  aa,ee,e21,e21s,ae,t2s, sigma, dd,dt2s,dt2s2 
real*8  pi, ax, by, cx 
integer  indx,mxindx 

pi=4.d0*datan(l .dO) 
c 

c  max;  number  of  iterations 
c 

mxindx=35 

c 

c  convergence  tolerance  on  consecutive  iterates 
c 

c  tol=l.e-18 

tol=l . e-6 
c 

c  convergence  tolerance  on  closeness  of  function  to  zero 
c 

c  tolf=l.e-18 

tolf=l . e-6 
c 

c  coefficients  of  transformation  matrix 
c 
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all*. 05 
al2*.05 
al3= . 05 
a21=.05 
a22=.05 
a23= . 05 

c  sigma  half  angle  of  shadow  cone 
sigma* . 26412*pi/180 . 
t2s=dtan(sigma)**2 

c  dd  distance  from  center  of  Earth  to  apex  of  shadow  cone 
dd=l .3836*10**6 
dt2s=dd*t2s 
c  ee  eccentricity 
ee=.001 

c  aa  semi  major  axis 
aa=10000000 
e21=l-ee*ee 
e21s=dsqrt  (e21) 
ae=aa*ee 

c  print  aa  ee  ',aa,ee 

c  print  *,'  e21  e21s  ae  ' ,e21 ,e21s,ae 

c 

c  initial  guess 

c 

xn=0 . 
yn=aa*e21 

print  *,'  xn  yn  initially  ' ,xn,yn 
c 

c  initialize  counter  of  iteration 

c 

indx=0 

c 

c  compute  coefficients  of  equation  to  solve 
c 

a=al2*al2+al3*al3-all*all*t2s 

b=a22*a22+a23*a23-a21*a21*t2s 


17 


c=2 . * (al2*a22+al3*a23-all*a21*t2s) 

dt2s2=2*dt2s 

d=dt2s2*all 

e=a21*dt2s2 

f=-d*dt2s 

print  a  b  c  d  e  f  ’ ,a,b,c,d,e,f 
5  continue 

call  fcn(a,b,c,d,e,f ,ae,e21,ax,by,cx,xn,yn,ypn,fn,fpn) 
print  indx  ypn  fn  fpn  ' ,indx,ypn,fn,fpn 
c 

c  is  function  close  to  zero  at  xn? 
c 

if (dabs(fn) .le.tolf)  go  to  20 
c 

c  compute  next  iterate 
c 

xnl=xn-fn/fpn 

ynl=e21s*dsqrt (aa*aa- (xnl+ae) **2) 
indx=indx+l 

print  *,indx,xn,xnl,yn,ynl,fn 
c 

c  check  for  closeness  of  iterates 
c 

if (dabs(xnl-xn) .le.tol)  go  to  10 
xn=xnl 
yn=ynl 
c 

c  check  if  max  number  of  iterates  exceeded 
c 

if  (indx. ge.mx indx)  go  to  30 
go  to  5 

10  print  convergence  iterates  close' ,xn,xnl 

go  to  40 

20  print  convergence  function  close  to  zero', xn,fn 

go  to  40 

30  print  no  convergence  -  max  number  of  iterates’ 
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40 


stop 
end 

subroutine  fcn(a,b .c^ejf  ,ae ,e21,ax,by }cx ,x,y ,yp,fn,fpn) 
real*8  x,y ,yp,fn,fpn 
real*8  a,b,c,d,e,f ,ax,by ,cx,ae,e21 
c 

c  evaluate  the  function 
c 

ax=a*x 

by=b*y 

cx=c*x 

fn=ax*x+by*y+cx*y+d*x+e*y+f 

c 

c  evaluates  the  first  derivative 
c 

yp=-e21* (x+ae) / y 

f pn=2 . *ax+2 . *by*yp+c*y+cx*yp+d+e*yp 

return 

end 
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